Importer les données

library(datasets)
data(iris)


Exploration de données

Avant de débuter des analyses, il est important de se familisariser avec son jeu de données afin d’avoir une idée de sa structure. Cette étape permet d’identifier des motifs, des tendances ou des relations préalablement aux tests statistiques, puis de vérifier la qualité des données. Elle permet aussi d’évaluer si les données ont besoin d’être transformées avant de procéder aux analyses.

Quelques fonctions de base

Pour vérifier le nombre de lignes et de colonnes de votre dataset:

dim(iris)
## [1] 150   5

Pour voir les 6 premières lignes du jeu de données:

head(iris)
##   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1          5.1         3.5          1.4         0.2  setosa
## 2          4.9         3.0          1.4         0.2  setosa
## 3          4.7         3.2          1.3         0.2  setosa
## 4          4.6         3.1          1.5         0.2  setosa
## 5          5.0         3.6          1.4         0.2  setosa
## 6          5.4         3.9          1.7         0.4  setosa

Pour connaître le nom des colonnes:

names(iris)
## [1] "Sepal.Length" "Sepal.Width"  "Petal.Length" "Petal.Width"  "Species"

Pour vérifier la classe de chaque variable dans le dataset:

str(iris)
## 'data.frame':    150 obs. of  5 variables:
##  $ Sepal.Length: num  5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
##  $ Sepal.Width : num  3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
##  $ Petal.Length: num  1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
##  $ Petal.Width : num  0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
##  $ Species     : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...

Pour obtenir certaines statistiques descriptives de base, telles que le minimum, le maximum puis la moyenne (variables continues), ainsi que le nombre d’observations (variables catégoriques):

summary(iris)
##   Sepal.Length    Sepal.Width     Petal.Length    Petal.Width   
##  Min.   :4.300   Min.   :2.000   Min.   :1.000   Min.   :0.100  
##  1st Qu.:5.100   1st Qu.:2.800   1st Qu.:1.600   1st Qu.:0.300  
##  Median :5.800   Median :3.000   Median :4.350   Median :1.300  
##  Mean   :5.843   Mean   :3.057   Mean   :3.758   Mean   :1.199  
##  3rd Qu.:6.400   3rd Qu.:3.300   3rd Qu.:5.100   3rd Qu.:1.800  
##  Max.   :7.900   Max.   :4.400   Max.   :6.900   Max.   :2.500  
##        Species  
##  setosa    :50  
##  versicolor:50  
##  virginica :50  
##                 
##                 
##

Pour oconnaître les niveaux d’une variable catégorique:

levels(iris$Species)
## [1] "setosa"     "versicolor" "virginica"


Visualisation d’un jeu de données

La fonction plot() est une fonction versatile qui permet de créer une grande variété de figures avec les données brutes. Dépendemment du type de donnée et des arguments fournis à la fonction, plot() peut générer plusieurs types de figures comme:

Une seule matrice de scatterplots figure peut être créée avec toutes les variables continues.

plot(iris)

Pour créer un scatterplot de la relation entre une variable spécifique et une autre, il suffit d’entrer celles-ci dans la fonction, séparées par le signe ~. Par exemple, si on veut visualiser la relation entre la longueur des sépales et celle des pétales dans le jeu de données iris, on peut l’écrire ainsi:

plot(iris$Sepal.Length~iris$Petal.Length)

À l’aide de la fonction boxplot(), on peut visualiser la dispersion d’une variable selon dles groupes d’une variable catégorique. Ce type de figure permet aussi d’identifier rapidement des valeurs aberrantes ou des anomalies.

boxplot(iris$Petal.Length~iris$Species)

La fonction hist() génère un histogramme de la variable qui y est précisée.

hist(iris$Sepal.Width)


Conseils pour des figures et tableaux réussis

Les figures sont des représentations visuelles des résultats. Elles rendent la lectures des résultats principaux plus facile et permettent de mettre en évidence des tendances ou motifs intéressants. On veut pouvoir la comprendre sans avoir à osciller entre la figure et le texte. Elles doivent contenir:

Les tableaux sont un bon choix pour présenter de l’information numérique détaillée. Ils présentent habituellement des résultats plus complexes qui seraient trop encombrants à inclure dans une figure ou dans le texte. Généralement, si les données ne peuvent être présentées en une ou deux phrases, un tableau est nécessaire. Les lignes et les colonnes doivent contenir le nom de la variable ainsi que l’unité de mesure. Un tableau résumé des statistiques peut inclure, par exemple, la moyenne, l’écart-type, les intervalles de confiances, les degrés de liberté, la valeur p et autres statistiques (comme la F value).

Légende d’une figure ou d’un tableau

Les légendes servent à compléter l’information qui est présentée dans le tableau ou la figure. On y retrouve entre autre les tailles d’échantillons (n), valeur p, des descriptions des abréviations, la méthode de collection, le nombre de réplicats, etc.

Annexe ou pas?

Les données ou les figures qui ne contribuent pas directement à l’histoire principale de votre rapport peuvent être rassemblées en annexes. Dans cette section, on peut retrouver, par exemple, des cartes du site d’échantillonnage, du matériel utilisé pour récolter les données, des tableaux avec davantage d’informations sur les modèles (pas de capture d’écran du summary!), etc.


Coder des graphiques avec ggplot (Sabrina)

Observer la stucture des données univariées

summary(iris)
##   Sepal.Length    Sepal.Width     Petal.Length    Petal.Width   
##  Min.   :4.300   Min.   :2.000   Min.   :1.000   Min.   :0.100  
##  1st Qu.:5.100   1st Qu.:2.800   1st Qu.:1.600   1st Qu.:0.300  
##  Median :5.800   Median :3.000   Median :4.350   Median :1.300  
##  Mean   :5.843   Mean   :3.057   Mean   :3.758   Mean   :1.199  
##  3rd Qu.:6.400   3rd Qu.:3.300   3rd Qu.:5.100   3rd Qu.:1.800  
##  Max.   :7.900   Max.   :4.400   Max.   :6.900   Max.   :2.500  
##        Species  
##  setosa    :50  
##  versicolor:50  
##  virginica :50  
##                 
##                 
##

ANOVA (Sabrina)

La fonction de base plot() peut être utile pour visualiser rapidement un jeu de données, par exemple avec un histogramme. Par contre, vous verrez que son utilisation peut devenir limitée lorsqu’il s’agit de réaliser des figures plus complexes ou simplement modifier certains paramètres graphiques.

C’est pour cela que nous suggérons d’utiliser la fonction ggplot(), du package “ggplot2”. Cette fonction permet de réaliser des graphiques de façon plus intuitive et permet de les mettre en page plus facilement qu’avec la fonction plot(). Elle est aussi mieux documentée, il est donc plus facile de comprendre et utiliser les différents aspects de la fonction. Google est d’ailleurs votre meilleur allié pour la réalisation de vos figures!

Contrairement à la fonction plot(), ggplot() fonctionne par couches. Une figure ggplot() commence avec la fonction ggplot(). Elle sert à “préparer” la figure: on spécifie le jeu de données à utiliser, puis on choisit les variables qui formeront nos axes.

La fonction ggplot() nécessite deux arguments: le dataset (jeu de données), puis l’argument aes(). Ce dernier nous permet d’assigner des variables du dataset aux composantes du graphique (par exemple, les axes x et y). Voici un exemple, toujours avec le jeu de données Iris. Nous allons tester si la largeur des sépales diffère entre les espèces.

(Revoir l’atelier 3 pour l’explication de l’ANOVA, incluant les postulats et critères d’utilisation, que je passe ici)

library("ggplot2")

head(iris)
##   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1          5.1         3.5          1.4         0.2  setosa
## 2          4.9         3.0          1.4         0.2  setosa
## 3          4.7         3.2          1.3         0.2  setosa
## 4          4.6         3.1          1.5         0.2  setosa
## 5          5.0         3.6          1.4         0.2  setosa
## 6          5.4         3.9          1.7         0.4  setosa
boxplot(Sepal.Width ~ Species,
  data = iris)

On produit d’abord le modèle linéaire, et le test d’ANOVA, avec Sepal.Width comme variable dépendante et Species comme variable indépendante.

modele.SW<-lm(Sepal.Width~Species,data=iris)

#Test anova:
anova(modele.SW)
## Analysis of Variance Table
## 
## Response: Sepal.Width
##            Df Sum Sq Mean Sq F value    Pr(>F)    
## Species     2 11.345  5.6725   49.16 < 2.2e-16 ***
## Residuals 147 16.962  0.1154                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(modele.SW)
## 
## Call:
## lm(formula = Sepal.Width ~ Species, data = iris)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -1.128 -0.228  0.026  0.226  0.972 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        3.42800    0.04804  71.359  < 2e-16 ***
## Speciesversicolor -0.65800    0.06794  -9.685  < 2e-16 ***
## Speciesvirginica  -0.45400    0.06794  -6.683 4.54e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3397 on 147 degrees of freedom
## Multiple R-squared:  0.4008, Adjusted R-squared:  0.3926 
## F-statistic: 49.16 on 2 and 147 DF,  p-value: < 2.2e-16
# [[[ Commentaire anova ]]]

Ensuite, on effectue le test post-hoc Tukey:

compSW<-aov(Sepal.Width~Species,data=iris)
TukeyHSD(compSW)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Sepal.Width ~ Species, data = iris)
## 
## $Species
##                        diff         lwr        upr     p adj
## versicolor-setosa    -0.658 -0.81885528 -0.4971447 0.0000000
## virginica-setosa     -0.454 -0.61485528 -0.2931447 0.0000000
## virginica-versicolor  0.204  0.04314472  0.3648553 0.0087802
summary(compSW)
##              Df Sum Sq Mean Sq F value Pr(>F)    
## Species       2  11.35   5.672   49.16 <2e-16 ***
## Residuals   147  16.96   0.115                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# [[[ Commentaire Tukey ]]]

On peut ensuite construire notre figure. D’abord, la fonction ggplot(). On spécifie le jeu de données, puis l’argument aes(), ici les axes x et y:

ggplot(iris, aes(x=Species, y=Sepal.Width))

Où sont les données? Il faut les ajouter!

C’est comme ça que ggplot() fonctionne. On crée la base de notre figure, puis on y ajoute les données à l’aide des fonctions geom_*. Par exemple, geom_points nous permet d’ajouter des données sous forme de scatterplots (donc des points), tandis que geom_line nous permet d’ajouter une ligne. Dans notre cas, nous voulons montrer nos données sous forme de boxplot, donc nous utilisons geom_boxplot. Il ne faut pas oublier d’ajouter un + après chaque fonction:

ggplot(iris, aes(x=Species, y=Sepal.Width)) +
geom_boxplot()

Pour ajouter/modifier d’autres éléments, il faut ajouter des couches. Ainsi, on peut changer les titres d’axes (labs pour labels), changer le thème de notre graphique (theme_bw est plus minimaliste), retirer la grille et centrer le titre:

ggplot(iris, aes(x=Species, y=Sepal.Width)) +
geom_boxplot()+
labs(   x="Espèce", 
    y="Largeur des sépales (cm)",
    title="Largeur des sépales chez les espèces d'iris")+
theme_bw()+
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

On peut aussi décider de changer les couleurs. Ce n’est pas toujours nécessaire ou pertinent, si ça n’ajoute pas d’information à la figure. Il faut faire attention de ne pas surcharger la figure avec le design esthétique, l’important c’est de focuser sur le message que l’on veut transmettre avec notre figure, et ne pas ajouter d’éléments qui sont distrayants. Mais pour l’exercice et se familiariser avec les paramètres graphiques, essayons:

ggplot(iris, aes(x=Species, y=Sepal.Width)) +
geom_boxplot(fill = "salmon", color = "red")+
labs(   x="Espèce", 
    y="Largeur des sépales (mm)",
    title="Largeur des sépales chez les espèces d'iris")+
theme_bw()+
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

On a donc changé le fill (intérieur) et color (trait) de nos boxplot, à même la fonction boxplot(). On pourrait aussi changer la couleur des boîtes en fonction de l’espèce d’iris (qui est aussi notre variable x). Dans ce cas, il faut ajouter un argument aes() à la fonction boxplot(), pour assigner notre variable à un aesthetic. Ce sera notre variable Species:

#Fill:

ggplot(iris, aes(x=Species, y=Sepal.Width)) +
geom_boxplot(aes(fill=Species), show.legend = FALSE)+
labs(   x="Espèce", 
    y="Largeur des sépales (mm)",
    title="Largeur des sépales chez les espèces d'iris")+
theme_bw()+
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) 

#Color:
ggplot(iris, aes(x=Species, y=Sepal.Width)) +
geom_boxplot(aes(color=Species), show.legend = FALSE)+
labs(   x="Espèce", 
    y="Largeur des sépales (mm)",
    title="Largeur des sépales chez les espèces d'iris")+
theme_bw()+
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

Il y a aussi des palettes de couleurs prédéfinies dans R, qu’on peut utiliser avec la fonction scale_fill_brewer (j’ai choisi la palette “Accent”):

ggplot(iris, aes(x=Species, y=Sepal.Width)) +
geom_boxplot(aes(fill=Species), show.legend=FALSE)+
scale_fill_brewer(palette="Accent")+
labs(   x="Espèce", 
    y="Largeur des sépales (mm)",
    title="Largeur des sépales chez les espèces d'iris")+
theme_bw()+
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

#Qu'est-ce qui arrive si on change Fill pour Color, dans boxplot()? Il faut aussi changer scale_fill_brewer() pour scale_color_brewer():

ggplot(iris, aes(x=Species, y=Sepal.Width)) +
geom_boxplot(aes(color=Species), show.legend=FALSE)+
scale_color_brewer(palette="Accent")+
labs(   x="Espèce", 
    y="Largeur des sépales (mm)",
    title="Largeur des sépales chez les espèces d'iris")+
theme_bw()+
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

Puisque ggplot() fonctionne par couche, on peut aussi choisir d’ajouter les données brutes à notre figure. Par exemple, on ajoute avant geom_boxplot() une autre couche de données avec geom_point(), dans lequel on peut spécifier, encore une fois, l’argument aes() pour assigner la variable qui sera colorée. L’argument position définit le niveau de “jitter”, c’est-à-dire des points un peu “éparpillés” sur l’axe x. Ensuite, on pourrait ajouter l’argument alpha=0.5 dans la fonction boxplot(), pour les rendre semi-transparents:

ggplot(iris, aes(x=Species, y=Sepal.Width)) +
geom_point(aes(color=Species), position = position_jitter(0.2), shape = 16, size = 2.5, show.legend=FALSE)+
scale_color_brewer(palette="Accent")+
geom_boxplot(alpha=0.5)+
labs(   x="Espèce", 
    y="Largeur des sépales (mm)",
    title="Largeur des sépales chez les espèces d'iris")+
theme_bw()+
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(hjust = 0.5))

On peut aussi ajouter des lettres pour montrer nos résultats de comparaison (voir https://statdoe.com/one-way-anova-and-box-plot-in-r/ pour la méthode). En gros, la fonction multcompLetters4 assigne les lettres en fonction des groupes du test de Tukey, puis on constuit une table qui contient ces lettres (LettresSW dans l’exemple), la variable Species, ainsi que la position de la lettre (en haut du quartile pour chaque boîte/valeur de Species).

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(multcompView)

compSW<-aov(Sepal.Width~Species,data=iris)
tukeySW <- TukeyHSD(compSW)
print(tukeySW)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Sepal.Width ~ Species, data = iris)
## 
## $Species
##                        diff         lwr        upr     p adj
## versicolor-setosa    -0.658 -0.81885528 -0.4971447 0.0000000
## virginica-setosa     -0.454 -0.61485528 -0.2931447 0.0000000
## virginica-versicolor  0.204  0.04314472  0.3648553 0.0087802
cld <- multcompLetters4(compSW, tukeySW)

LettresSW <- group_by(iris, Species) %>%
  summarise(mean=mean(Sepal.Width), quant = quantile(Sepal.Width, probs = 0.75)) %>%
  arrange(desc(mean))

cld <- as.data.frame.list(cld$Species)
LettresSW$cld <- cld$Letters

print(LettresSW)
## # A tibble: 3 × 4
##   Species     mean quant cld  
##   <fct>      <dbl> <dbl> <chr>
## 1 setosa      3.43  3.68 a    
## 2 virginica   2.97  3.18 b    
## 3 versicolor  2.77  3    c

On ajoute nos lettres au boxplot avec la fonction geom_text(), avec comme argument aes() la table que l’on vient de créer. On peut aussi ajouter la statistique de test avec la fonction annotate(), en utilisant les coordonnées x et y dans notre figure.

ggplot(iris, aes(x=Species, y=Sepal.Width)) +
geom_point(aes(color=Species), position = position_jitter(0.1), shape = 16, size = 2.5, show.legend=FALSE)+
scale_color_brewer(palette="Accent")+
geom_boxplot(alpha=0.5)+
labs(   x="Espèce", 
    y="Largeur des sépales (mm)",
    title="Largeur des sépales chez les espèces d'iris")+
theme_bw()+
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(hjust = 0.5))+
geom_text(data = LettresSW, aes(x = Species, y = quant, label = cld), size = 5, vjust=-1, hjust =-5)+
annotate("text", x = 3, y = 4.5, label = "ANOVA, F(2)=49.16, p<0.001")

Pour bien faire, on pourrait même changer les étiquettes en x. Avec la fonction scale_x_discrete, j’ajoute “I.” pour iris à chaque étiquette et je les mets en italique, puisque ce sont des noms latins, avec l’argument axis.text.x dans la fonction theme().

ggplot(iris, aes(x=Species, y=Sepal.Width)) +
geom_point(aes(color=Species), position = position_jitter(0.1), shape = 16, size = 2.5, show.legend=FALSE)+
scale_color_brewer(palette="Accent")+
geom_boxplot(alpha=0.5)+
labs(   x="Espèce", 
    y="Largeur des sépales (mm)",
    title="Largeur des sépales chez les espèces d'iris")+
theme_bw()+
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(hjust = 0.5), axis.text.x = element_text(face="italic"))+
geom_text(data = LettresSW, aes(x = Species, y = quant, label = cld), size = 5, vjust=-1, hjust =-5)+
annotate("text", x = 3, y = 4.5, label = "ANOVA, F(2)=49.16, p<0.001")+
scale_x_discrete(labels=c("I. setosa", "I. versicolor", "I. virginica"))

Régression linéaire (Sabrina)

Autres modèles univariés… (Sabrina, à voir si c’est pertinent)

Ordinations

Comme pour les analyses univariées, cet atelier est seulement axé sur la création de graphiques. Pour une introduction aux analyses multivariées (ordinations, PERMANOVA, Procrustes), consulter l’atelier 6.

Importer les données multivariées

library(vegan)
## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.6-4
data(dune) # matrice de communautés
data(dune.env) # données environnementales

Ce jeu de données contient une matrice d’abondance (dune) d’espèces végétales avec leur classe de couverture pour 20 sites. Le dataframe (dune.env) contient 5 variables denvironnementales.

PCA

Préparation des données

Transformation Hellinger

dune.hel <- decostand(dune, method="hellinger")

Faire une PCA (avec la fonction prcomp cette fois-ci, car l’objet qu’on va créer peut être “lu” plus tard par la fonction ggord)

pca.dune <- prcomp(dune.hel)
summary(pca.dune)
## Importance of components:
##                           PC1    PC2    PC3    PC4     PC5     PC6     PC7
## Standard deviation     0.4122 0.3333 0.2374 0.2097 0.20242 0.17376 0.14506
## Proportion of Variance 0.3057 0.1998 0.1013 0.0791 0.07371 0.05431 0.03785
## Cumulative Proportion  0.3057 0.5054 0.6068 0.6859 0.75960 0.81391 0.85177
##                            PC8     PC9    PC10    PC11    PC12    PC13    PC14
## Standard deviation     0.13392 0.11745 0.10441 0.09995 0.09550 0.07985 0.07188
## Proportion of Variance 0.03226 0.02481 0.01961 0.01797 0.01641 0.01147 0.00929
## Cumulative Proportion  0.88403 0.90885 0.92846 0.94643 0.96283 0.97430 0.98360
##                           PC15    PC16    PC17    PC18    PC19      PC20
## Standard deviation     0.05424 0.05240 0.04063 0.03373 0.02536 1.816e-17
## Proportion of Variance 0.00529 0.00494 0.00297 0.00205 0.00116 0.000e+00
## Cumulative Proportion  0.98889 0.99383 0.99680 0.99884 1.00000 1.000e+00

Dans l’atelier 6, nous avons fait nos PCA avec la fonction rda. Lorsqu’on compare le sommaire des deux analyses, on voit que les deux fonctions donnent les mêmes résultats.

pca2.dune <- rda(dune.hel)
summary(pca2.dune)
## 
## Call:
## rda(X = dune.hel) 
## 
## Partitioning of variance:
##               Inertia Proportion
## Total          0.5559          1
## Unconstrained  0.5559          1
## 
## Eigenvalues, and their contribution to the variance 
## 
## Importance of components:
##                          PC1    PC2     PC3     PC4     PC5     PC6     PC7
## Eigenvalue            0.1699 0.1111 0.05634 0.04397 0.04097 0.03019 0.02104
## Proportion Explained  0.3057 0.1998 0.10135 0.07910 0.07371 0.05431 0.03785
## Cumulative Proportion 0.3057 0.5054 0.60679 0.68590 0.75960 0.81391 0.85177
##                           PC8     PC9    PC10     PC11    PC12     PC13
## Eigenvalue            0.01794 0.01379 0.01090 0.009991 0.00912 0.006376
## Proportion Explained  0.03226 0.02481 0.01961 0.017972 0.01641 0.011469
## Cumulative Proportion 0.88403 0.90885 0.92846 0.946428 0.96283 0.974303
##                           PC14     PC15     PC16     PC17     PC18      PC19
## Eigenvalue            0.005166 0.002942 0.002745 0.001651 0.001138 0.0006432
## Proportion Explained  0.009293 0.005292 0.004939 0.002969 0.002047 0.0011570
## Cumulative Proportion 0.983597 0.988888 0.993827 0.996796 0.998843 1.0000000
## 
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
## * General scaling constant of scores:  1.80276 
## 
## 
## Species scores
## 
##                 PC1       PC2       PC3       PC4        PC5        PC6
## Achimill -0.2241355  0.057048 -0.068674 -0.173281 -0.0498964 -0.0006022
## Agrostol  0.4142446 -0.180068  0.020078 -0.020756  0.0029054  0.0193453
## Airaprae -0.0350694  0.154644  0.116317 -0.039999 -0.1384374 -0.0376930
## Alopgeni  0.1386550 -0.318444  0.186419  0.001329  0.0030171 -0.0215368
## Anthodor -0.1962910  0.238579  0.095768 -0.179434 -0.0380470 -0.0471885
## Bellpere -0.1148332 -0.054380 -0.057640  0.064957  0.0010453  0.0896219
## Bromhord -0.1410948 -0.053134 -0.049586 -0.057365  0.0210053  0.0938974
## Chenalbu  0.0121539 -0.025919  0.036102 -0.015353 -0.0038179  0.0173353
## Cirsarve -0.0026599 -0.031975  0.004872  0.021813 -0.0182940  0.0064877
## Comapalu  0.1100812  0.048635 -0.061981 -0.028257  0.0102828  0.1064664
## Eleopalu  0.3796078  0.069099 -0.185862 -0.064110  0.0144632  0.0035959
## Elymrepe -0.1456860 -0.248314 -0.094249  0.010195 -0.1482516 -0.0563573
## Empenigr  0.0004895  0.057213  0.063829  0.026314 -0.0393467 -0.0066958
## Hyporadi -0.0559119  0.197829  0.139068  0.036814 -0.1427093 -0.0326104
## Juncarti  0.2478729 -0.012497 -0.092252  0.005525  0.0265396 -0.1569375
## Juncbufo  0.0190348 -0.121457  0.179346 -0.051491  0.0931930 -0.0211675
## Lolipere -0.3366259 -0.173707 -0.206664  0.154192  0.0092155 -0.0325721
## Planlanc -0.2547353  0.197554 -0.034482 -0.039980  0.1469090 -0.0256698
## Poaprat  -0.2808972 -0.142814 -0.082463  0.084918 -0.0502270 -0.0278725
## Poatriv  -0.1317430 -0.363370  0.059442 -0.137436  0.0581366 -0.0136469
## Ranuflam  0.2783344  0.024399 -0.081013 -0.052662  0.0009058  0.0274683
## Rumeacet -0.1079283 -0.024094  0.054027 -0.122507  0.2108960 -0.0801798
## Sagiproc  0.0409392 -0.082141  0.246528  0.127747 -0.0148906 -0.0086761
## Salirepe  0.0651846  0.156446  0.007210  0.136628 -0.0166452 -0.0581183
## Scorautu -0.0562239  0.158686  0.094988  0.095244  0.0265794  0.1029567
## Trifprat -0.1001443  0.028383 -0.024352 -0.095019  0.1389927 -0.0417313
## Trifrepe -0.0471194 -0.006999  0.052812  0.045812  0.1395303  0.2537096
## Vicilath -0.0541878  0.052774 -0.023104  0.114929  0.0381791  0.0343033
## Bracruta  0.0852752  0.107955 -0.003783  0.181534  0.2198174 -0.1398609
## Callcusp  0.2045115  0.057120 -0.104072 -0.065216 -0.0146101  0.0579472
## 
## 
## Site scores (weighted sums of species scores)
## 
##          PC1      PC2        PC3      PC4       PC5       PC6
## 1  -0.423423 -0.42754 -6.533e-01  0.01583 -0.761705 -0.463469
## 2  -0.356333 -0.36816 -1.991e-01 -0.08370 -0.355682  0.630146
## 3  -0.080900 -0.55215 -5.753e-02  0.27275 -0.217514  0.016724
## 4  -0.041005 -0.49292  7.510e-02  0.33626 -0.282018  0.100013
## 5  -0.458003  0.05220 -1.432e-01 -0.54110  0.351891 -0.117334
## 6  -0.389342  0.20337 -1.022e-01 -0.30549  0.758789 -0.245696
## 7  -0.451813  0.06864 -6.835e-02 -0.41821  0.585536 -0.138732
## 8   0.295113 -0.29375 -6.362e-02  0.19828  0.010903 -0.185034
## 9   0.035973 -0.49264  2.714e-01  0.07945  0.086166 -0.419257
## 10 -0.453094  0.13639 -2.391e-01 -0.12201  0.144355  0.462946
## 11 -0.273214  0.29633  3.931e-05  0.87663  0.126762  0.097096
## 12  0.246622 -0.33173  9.205e-01 -0.03527  0.493526 -0.017653
## 13  0.226908 -0.48390  6.740e-01 -0.28663 -0.071279  0.323642
## 14  0.570818  0.25288 -3.165e-01 -0.32869 -0.085439  1.196310
## 15  0.654415  0.28845 -3.733e-01  0.01035  0.196968  0.002255
## 16  0.720030 -0.10689 -2.727e-01 -0.35992  0.007984 -0.504518
## 17 -0.317465  0.75272  3.395e-01 -0.64284 -0.803236 -0.262492
## 18 -0.201131  0.39818 -2.007e-01  0.89875  0.365682  0.086314
## 19  0.006263  0.73205  8.167e-01  0.33669 -0.503443 -0.085674
## 20  0.689579  0.36848 -4.077e-01  0.09888 -0.048247 -0.475587

PERMANOVA. On va tester si la composition végétale diffère entre différents type d’aménagement (facteur; variable quantitative) des sites.

dune.dist <- vegdist(decostand(dune, method='hellinger'), method='euclid')
disper.dune <- betadisper(dune.dist, dune.env$Management)
anova(disper.dune)
## Analysis of Variance Table
## 
## Response: Distances
##           Df  Sum Sq  Mean Sq F value Pr(>F)
## Groups     3 0.17812 0.059373  1.5729 0.2349
## Residuals 16 0.60397 0.037748
permanova.dune <- adonis2(dune.dist ~ Management, data = dune.env)
permanova.dune
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = dune.dist ~ Management, data = dune.env)
##            Df SumOfSqs      R2      F Pr(>F)   
## Management  3   3.1672 0.29986 2.2842  0.005 **
## Residual   16   7.3950 0.70014                 
## Total      19  10.5621 1.00000                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Le test montre que le type d’aménagement a un effet statistiquement significatif sur la composition en espèces végétales des sites.

ggord + ggplot

On va créer un graphique pour présenter la PCA avec les différents types d’aménagement.

Extraire les scores des sites et des espèces, à partir de l’objet créer par prcomp. Cela va nous permettre de mettre éventuellement seulement les espèces désirées sur le graphique.

site.scores<-pca.dune$x 
site.scores.df <- data.frame(site.scores)
sp.scores<-pca.dune$rotation
sp.scores.df <- data.frame(sp.scores)

Ouverture des librairies nécessaires à la création du graphique. Si problème avec installation de la librairie ggord : https://rdrr.io/github/fawda123/ggord/

library(ggord)
library(ggplot2)
library(ggrepel) # pour éviter overlap du texte

On va d’abord commencer le graphique avec la fonctions ggord du package du même nom.

Le graphique créé par la fonction ggord est un objet ggplot. Cet objet (ici appelé pca.plot) peut être personnalisé à même la fonction ggord, et peut aussi servir de base sur laquelle ajouter des couches avec ggplot.

pca.plot <- ggord(pca.dune, xlims=c(-1,1), ylims=c(-1,1)) # xlims et ylims définissent les limites dans l'affichage des axe (les valeurs entre parenthèses correspondent à l'échelle numérique des axes). On ajuste au fur et à mesure, au besoin.
pca.plot # À la fin de chaque chunk, l'objet graphique est appelé pour qu'on visualise le rendu à chaque fois.

On remarque que la fonctionggord inclut par défaut la proportion de la variance expliquée par les axes. C’est super!

Par défaut, la fonction ggord va mettre en graphique les deux premiers axes de l’ordination. Mais si on s’intéresse à différents axes, on peut spécifier lesquels mettre en graphique.

pca.plot <- ggord(pca.dune, xlims=c(-1,1), ylims=c(-1,1),
                  axes=c(1,3)) # Ici, on met en graphique les axes 1 et 3.
pca.plot

On revient aux deux premiers axes. Faisons quelques petits changements. Les changements sont marqués d’un “#” à côté de la ligne.

pca.plot <- ggord(pca.dune, xlims=c(-1.2,1.4), ylims=c(-0.75,0.9), # ajustement des limites
                  axes=c(1,2), # retour aux axes 1 et 2
                  size=3, # diminution de la taille des points (sites)
                  arrow=NULL,# enlever les flèches (qui viennent par défaut) pour les espèces
                  labcol="forestgreen") # changer la couleur du texte pour les espèces
pca.plot

Faisons d’autres changements!

pca.plot <- ggord(pca.dune, xlims=c(-1.2,1.4), ylims=c(-0.75,0.9), 
                  axes=c(1,2), 
                  size=3, 
                  obslab=TRUE, # pour passer de points pour les sites (par défaut) à texte (nom des sites)
                  arrow=NULL,
                  labcol="forestgreen") # changer la couleur du texte pour les espèces
pca.plot

Il y a beaucoup d’espèces affichées au milieu et c’est illisible. En plus, les espèces au milieu ne contribuent pas beaucoup à la distinction compositionnelle des sites. On va donc garder seulement les quelques espèces dont les scores ont la valeur absolue la plus élevée, pour chacun des deux axes, qu’on met en graphique (ici PC1 et PC2).

D’abord, on va créer un nouveau dataframe contenant seulement ces espèces, en utilisant quelques fonctions de la librairie dplyr.

A <- top_n(sp.scores.df, 3, PC1) # sélectionne les 3 espèces dont le score est le plus élevé sur l'axe 1.
B <- top_n(sp.scores.df, -3, PC1) # sélectionne les 3 espèces dont le score est le plus bas sur l'axe 1.
C <- top_n(sp.scores.df, 3, PC2) # sélectionne les 3 espèces dont le score est le plus élevé sur l'axe 2.
D <- top_n(sp.scores.df, -3, PC2) # sélectionne les 3 espèces dont le score est le plus bas sur l'axe 2.
sp.scores.skim <- bind_rows(A,B,C,D) # merge les sub dataframe A, B, C et D.
sp.scores.skim <- distinct(sp.scores.skim) # on ne garde que les rangées (espèces) qui sont uniques (pour éviter qu'une même espèce s'y retrouve en copies)

Mise en graphique de ces espèces seulement

pca.plot <- ggord(pca.dune, xlims=c(-1.2,1.4), ylims=c(-0.75,0.9), 
                  axes=c(1,2), 
                  size=3, 
                  obslab=TRUE,
                  arrow=NULL,
                  txt=NULL, # enlever le texte (par défaut) pour les espèces
                  addpts=sp.scores.skim, # Ajouter les espèces aux plus hauts et faibles scores
                  addsize=3, # taille des points des espèces
                  addcol="forestgreen") # couleur des espèces
                  # on enlève l'argument labcol
pca.plot

pca.plot <- ggord(pca.dune, xlims=c(-1.2,1.4), ylims=c(-0.75,0.9), 
                  axes=c(1,2), 
                  size=3, 
                  obslab=TRUE,
                  arrow=NULL,
                  txt=NULL, 
                  addpts=sp.scores.skim, 
                  addsize=4, # augmenter la taille du texte des espèces
                  ptslab=TRUE, # changer de points à texte pour les espèces
                  addcol="forestgreen")
pca.plot

On va ajouter les résultats d’un envfit. D’abord, on fait le test.

(pca.envfit <- envfit(pca.dune, dune.env))
## 
## ***VECTORS
## 
##        PC1     PC2     r2 Pr(>r)   
## A1 0.97222 0.23409 0.3622  0.008 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
## 
## ***FACTORS:
## 
## Centroids:
##                  PC1     PC2
## Moisture1    -0.3635  0.0257
## Moisture2    -0.2224 -0.0314
## Moisture4     0.1408 -0.3321
## Moisture5     0.4504  0.0872
## ManagementBF -0.3597  0.0173
## ManagementHF -0.1930 -0.0745
## ManagementNM  0.2330  0.3751
## ManagementSF  0.1077 -0.3217
## UseHayfield  -0.0994  0.2242
## UseHaypastu  -0.0203 -0.2180
## UsePasture    0.1716  0.0350
## Manure0       0.2330  0.3751
## Manure1      -0.2294 -0.0161
## Manure2      -0.2385 -0.0895
## Manure3       0.1969 -0.1644
## Manure4      -0.1812 -0.3955
## 
## Goodness of fit:
##                r2 Pr(>r)    
## Moisture   0.5366  0.001 ***
## Management 0.4614  0.001 ***
## Use        0.1794  0.131    
## Manure     0.4531  0.010 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999

Pour incorporer les vecteurs résultant de ce test, il faut d’abord les extraire et les mettre dans un dataframe.

vect <- data.frame(pca.envfit[["vectors"]][["arrows"]])
vect$variable.env<-rownames(vect)
pca.plot <- ggord(pca.dune, xlims=c(-1.2,1.4), ylims=c(-0.75,0.9), 
                  axes=c(1,2), 
                  size=3, 
                  obslab=TRUE,
                  arrow=NULL,
                  txt=NULL, 
                  addpts=sp.scores.skim, 
                  addsize=4, 
                  ptslab=TRUE, 
                  addcol="forestgreen")
pca.plot + 
  coord_fixed() + # ajout des vecteurs de envfit
  geom_segment(data = vect,
               aes(x = 0, xend = PC1, y = 0, yend = PC2),
               arrow = arrow(length = unit(0.5, "cm")), colour = "darkgrey") +
  geom_text_repel(data = vect, aes(x = PC1, y = PC2, label = variable.env, size = 3), show.legend = FALSE)
## Coordinate system already present. Adding new coordinate system, which will
## replace the existing one.

Maintenant, on va ajouter des ellipses autour des sites appartenant à des groupes similaires. On va aussi revenir à des symboles au lieu du nom des sites et changer la couleur du texte des espèces.

pca.plot <- ggord(pca.dune, xlims=c(-1.2,1.4), ylims=c(-0.75,0.9),
                  axes=c(1,2),
                  grp_in=dune.env$Management, # ajout des ellipses en fonction du type de Management
                  grp_title = "Management",
                  ellipse_pro=0.95, # valeur de confiance pour les ellipses
                  alpha_el=0.3, # transparence des ellipses
                  size=3, 
                  obslab=FALSE, # changement de texte à points pour sites
                  arrow=NULL,
                  txt=NULL, 
                  addpts=sp.scores.skim, 
                  addsize=4, 
                  ptslab=TRUE,
                  addcol="black", # couleur des espèces
                  alpha = 1)
pca.plot + 
    scale_shape_manual('Groups', values = c(15,16,17,18)) + # Symboles pour les sites en fonction de leur groupe. Assurer vous d'avoir un nombre de symboles correspondant au nombre de groupes.
    coord_fixed() +
    geom_segment(data = vect,
               aes(x = 0, xend = PC1, y = 0, yend = PC2),
               arrow = arrow(length = unit(0.5, "cm")), colour = "darkgrey") +
    geom_text_repel(data = vect, aes(x = PC1, y = PC2, label = variable.env, size = 3), show.legend = FALSE)
## Scale for shape is already present.
## Adding another scale for shape, which will replace the existing scale.
## Coordinate system already present. Adding new coordinate system, which will
## replace the existing one.

Voici une légende indiquant quels symboles sont associées aux différentes valeurs numériques de la couche scale_shape_manual.

Changer les couleurs des groupes On va extraire les codes des couleurs d’une palette de RColorBrewer avec la fonction brewer.pal

library(RColorBrewer)
pca.colors <- brewer.pal(n = 4, name = "Dark2") # Hexadecimal color specification

On va incorporer ces codes de couleurs dans un argument. On déplace la légende vers le haut.

pca.plot <- ggord(pca.dune, xlims=c(-1.2,1.4), ylims=c(-0.75,0.9),
                  axes=c(1,2),
                  grp_in=dune.env$Management, 
                  grp_title = "Management",
                  cols=pca.colors, # ajout des couleurs choisies
                  ellipse_pro=0.95,
                  alpha_el=0.3,
                  size=3, 
                  obslab=FALSE, 
                  arrow=NULL,
                  txt=NULL, 
                  addpts=sp.scores.skim, 
                  addsize=4, 
                  ptslab=TRUE,
                  addcol="black",
                  alpha = 1)
pca.plot <- pca.plot + # ici j'enregistre l'ensemble des arguments précédents et suivants dans mon objet pca.plot.
    scale_shape_manual('Groups', values = c(15,16,17,18)) + 
    coord_fixed() +
    geom_segment(data = vect,
               aes(x = 0, xend = PC1, y = 0, yend = PC2),
               arrow = arrow(length = unit(0.5, "cm")), colour = "darkgrey") +
    geom_text_repel(data = vect, aes(x = PC1, y = PC2, label = variable.env, size = 3), show.legend = FALSE)+
    theme(legend.position = 'top') # mettre la légende en haut
## Scale for shape is already present.
## Adding another scale for shape, which will replace the existing scale.
## Coordinate system already present. Adding new coordinate system, which will
## replace the existing one.
pca.plot

Créer un thème custom pour ggplot.

super_theme <- theme(
        panel.background = element_blank(),
        #panel.border = element_blank(),
        panel.grid = element_blank(),
        axis.line = element_line("gray25"),
        text = element_text(size = 12),
        axis.text = element_text(size = 10, colour = "gray25"),
        axis.title = element_text(size = 14, colour = "gray25"),
        legend.title = element_text(size = 14),
        legend.text = element_text(size = 14),
        legend.key = element_blank())

Ajouter ce thème au graphique.

pca.plot + super_theme

On va rajouter du texte!

pca.plot + annotate("text", x = 1.1, y = 0.7, label = "PERMANOVA\nR2=30%\n p=0.007") +
    super_theme

Satisfait? On sauvegarde l’image en png. Par défaut, ggplot enregistre le dernier graphique créé. Ici on copie-colle le code du graphique ci haut, puis sauvegarde ce graphique avec la fonction ggsave dans un dossier qu’on a créé (figures/).

pca.plot + annotate("text", x = 1.1, y = 0.7, label = "PERMANOVA\nR2=30%\n p=0.007") +
    super_theme

ggsave("pca.png", path="./figures/")
## Saving 7 x 5 in image

Cet atelier n’est pas fait pour vous montrer à utiliser une esthétique particulière, mais bien pour vous familiariser avec les options (infinies…) permettant de personnaliser un graphique d’ordination. Vous allez voir des figures d’ordination dans la littérature. Essayez de reproduire l’esthétique qui vous semble la plus appropriée… ou, tout simplement, amusez-vous les artistes!!

CA

Préparation des données

Faire une CA avec la fonction cca.

ca.dune <- cca(dune)
ca.summary <- summary(ca.dune)
summary(ca.dune)
## 
## Call:
## cca(X = dune) 
## 
## Partitioning of scaled Chi-square:
##               Inertia Proportion
## Total           2.115          1
## Unconstrained   2.115          1
## 
## Eigenvalues, and their contribution to the scaled Chi-square 
## 
## Importance of components:
##                          CA1    CA2    CA3     CA4     CA5     CA6     CA7
## Eigenvalue            0.5360 0.4001 0.2598 0.17598 0.14476 0.10791 0.09247
## Proportion Explained  0.2534 0.1892 0.1228 0.08319 0.06844 0.05102 0.04372
## Cumulative Proportion 0.2534 0.4426 0.5654 0.64858 0.71702 0.76804 0.81175
##                           CA8     CA9    CA10    CA11    CA12    CA13     CA14
## Eigenvalue            0.08091 0.07332 0.05630 0.04826 0.04125 0.03523 0.020529
## Proportion Explained  0.03825 0.03466 0.02661 0.02282 0.01950 0.01665 0.009705
## Cumulative Proportion 0.85000 0.88467 0.91128 0.93410 0.95360 0.97025 0.979955
##                           CA15     CA16     CA17     CA18     CA19
## Eigenvalue            0.014911 0.009074 0.007938 0.007002 0.003477
## Proportion Explained  0.007049 0.004290 0.003753 0.003310 0.001644
## Cumulative Proportion 0.987004 0.991293 0.995046 0.998356 1.000000
## 
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
## 
## 
## Species scores
## 
##               CA1      CA2      CA3       CA4       CA5      CA6
## Achimill -0.90859  0.08461 -0.58636 -0.008919 -0.660183 -0.18877
## Agrostol  0.93378 -0.20651  0.28165  0.024293 -0.139326 -0.02256
## Airaprae -1.00434  3.06749  1.33773  0.194305 -1.081813 -0.53699
## Alopgeni  0.40088 -0.61839  0.85013  0.346740  0.016509  0.10169
## Anthodor -0.96676  1.08361 -0.17188  0.459788 -0.607533 -0.30425
## Bellpere -0.50018 -0.35503 -0.15239 -0.704153 -0.058546  0.07308
## Bromhord -0.65762 -0.40634 -0.30685 -0.496751 -0.561358  0.07004
## Chenalbu  0.42445 -0.84402  1.59029  1.248755 -0.207480  0.87566
## Cirsarve -0.05647 -0.76398  0.91793 -1.175919 -0.384024 -0.13985
## Comapalu  1.91690  0.52150 -1.18215 -0.021738 -1.359988  1.31207
## Eleopalu  1.76383  0.34562 -0.57336 -0.002976 -0.332396 -0.14688
## Elymrepe -0.37074 -0.74148  0.26238 -0.566308 -0.270122 -0.72624
## Empenigr -0.69027  3.26420  1.95716 -0.176936 -0.073518 -0.16083
## Hyporadi -0.85408  2.52821  1.13951 -0.175115 -0.311874  0.11177
## Juncarti  1.27580  0.09963 -0.09320  0.005536  0.289410 -0.78247
## Juncbufo  0.08157 -0.68074  1.00545  1.078390  0.268360  0.24168
## Lolipere -0.50272 -0.35955 -0.21821 -0.474727  0.101494 -0.01594
## Planlanc -0.84058  0.24886 -0.78066  0.371149  0.271377  0.11989
## Poaprat  -0.38919 -0.32999 -0.02015 -0.358371  0.079296 -0.05165
## Poatriv  -0.18185 -0.53997  0.23388  0.178834 -0.155342 -0.07584
## Ranuflam  1.55886  0.30700 -0.29765  0.046974 -0.008747 -0.14744
## Rumeacet -0.65289 -0.25525 -0.59728  1.160164  0.255849 -0.32730
## Sagiproc  0.00364  0.01719  1.11570  0.066981  0.186654  0.32463
## Salirepe  0.61035  1.54868  0.04970 -0.607136  1.429729 -0.55183
## Scorautu -0.19566  0.38884  0.03975 -0.130392  0.141232  0.23717
## Trifprat -0.88116 -0.09792 -1.18172  1.282429  0.325706 -0.33388
## Trifrepe -0.07666 -0.02032 -0.20594  0.026462 -0.186748  0.53957
## Vicilath -0.61893  0.37140 -0.46057 -1.000375  1.162652  1.44971
## Bracruta  0.18222  0.26477 -0.16606  0.064009  0.576334  0.07741
## Callcusp  1.95199  0.56743 -0.85948 -0.098969 -0.556737  0.23282
## 
## 
## Site scores (weighted averages of species scores)
## 
##         CA1        CA2      CA3      CA4      CA5      CA6
## 1  -0.81167 -1.0826714 -0.14479 -2.10665 -0.39287 -1.83462
## 2  -0.63268 -0.6958357 -0.09708 -1.18695 -0.97686  0.06575
## 3  -0.10148 -0.9128732  0.68815 -0.68137 -0.08709 -0.28678
## 4  -0.05647 -0.7639784  0.91793 -1.17592 -0.38402 -0.13985
## 5  -0.95293 -0.1846015 -0.95609  0.86853 -0.34552 -0.98333
## 6  -0.85633 -0.0005408 -1.39735  1.59909  0.65494 -0.19386
## 7  -0.87149 -0.2547040 -0.86830  0.90468  0.17385 -0.03446
## 8   0.76268 -0.2968459  0.35648 -0.10772  0.17507 -0.36444
## 9   0.09693 -0.7864314  0.86492  0.40090  0.28704 -1.02783
## 10 -0.87885 -0.0353136 -0.82987 -0.68053 -0.75438  0.81070
## 11 -0.64223  0.4440332 -0.17371 -1.09684  1.37462  2.00626
## 12  0.28557 -0.6656161  1.64423  1.71496  0.65381  1.17376
## 13  0.42445 -0.8440195  1.59029  1.24876 -0.20748  0.87566
## 14  1.91996  0.5351062 -1.39863 -0.08575 -2.21317  2.43044
## 15  1.91384  0.5079036 -0.96567  0.04227 -0.50681  0.19370
## 16  2.00229  0.1090627 -0.33414  0.33760 -0.50097 -0.76159
## 17 -1.47545  2.7724102  0.40859  0.75117 -2.59425 -1.10122
## 18 -0.31241  0.6328355 -0.66501 -1.12728  2.65575  0.97565
## 19 -0.69027  3.2642026  1.95716 -0.17694 -0.07352 -0.16083
## 20  1.94438  1.0688809 -0.66595 -0.55317  1.59606 -1.70292

Extraction de la proportion de la variance expliquée par les deux premiers axes (pour s’y référer plus rapidement).

prop.expl.ca1 <- ca.summary$cont$importance[2, "CA1"]
prop.expl.ca2 <- ca.summary$cont$importance[2, "CA2"]

envfit

(ca.envfit <- envfit(ca.dune, dune.env))
## 
## ***VECTORS
## 
##         CA1      CA2     r2 Pr(>r)  
## A1 0.998160 0.060614 0.3104  0.044 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
## 
## ***FACTORS:
## 
## Centroids:
##                  CA1     CA2
## Moisture1    -0.7484 -0.1423
## Moisture2    -0.4652 -0.2156
## Moisture4     0.1827 -0.7315
## Moisture5     1.1143  0.5708
## ManagementBF -0.7258 -0.1413
## ManagementHF -0.3867 -0.2960
## ManagementNM  0.6517  1.4405
## ManagementSF  0.3376 -0.6761
## UseHayfield  -0.2861  0.6488
## UseHaypastu  -0.0735 -0.5602
## UsePasture    0.5163  0.0508
## Manure0       0.6517  1.4405
## Manure1      -0.4639 -0.1738
## Manure2      -0.5872 -0.3600
## Manure3       0.5187 -0.3172
## Manure4      -0.2059 -0.8775
## 
## Goodness of fit:
##                r2 Pr(>r)    
## Moisture   0.4113  0.005 ** 
## Management 0.4441  0.001 ***
## Use        0.1845  0.076 .  
## Manure     0.4552  0.006 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999

Pour incorporer les vecteurs résultant de ce test, il faut d’abord les extraire et les mettre dans un dataframe.

vect.ca.env <- data.frame(ca.envfit[["vectors"]][["arrows"]])
vect.ca.env$variable.env<-rownames(vect.ca.env)

PERMANOVA. On va tester si la composition végétale diffère entre différents type d’aménagement (facteur; variable quantitative) des sites.

dune.dist.chi2 <- vegdist(dune, method='chisq')
disper.dune.chi2 <- betadisper(dune.dist.chi2, dune.env$Management)
anova(disper.dune.chi2)
## Analysis of Variance Table
## 
## Response: Distances
##           Df Sum Sq Mean Sq F value   Pr(>F)   
## Groups     3 3.0277  1.0092  5.9124 0.006493 **
## Residuals 16 2.7311  0.1707                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
permanova.dune.ca <- adonis2(dune.dist.chi2 ~ Management, data = dune.env)
permanova.dune.ca
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = dune.dist.chi2 ~ Management, data = dune.env)
##            Df SumOfSqs      R2      F Pr(>F)   
## Management  3   12.399 0.25508 1.8263  0.006 **
## Residual   16   36.208 0.74492                 
## Total      19   48.606 1.00000                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Ici, le résultat du test de dispersion indique que la différence dans la composition végétale entre les types d’aménagement qui est “perçue” par la PERMANOVA pourrait être du à la différence dans la variance interne en composition des groupes.

Extraire les scores de CA avec la fonction scores

site.scores.dune.ca <- scores(ca.dune, display="sites", choices=c(1,2)) # ici on garde que les scores des axes 1 et 2 mais on aurait pu extraire tous les scores aussi.
site.scores.dune.ca <- data.frame(site.scores.dune.ca)
Management <- dune.env$Management
site.scores.ca <- cbind(site.scores.dune.ca, Management)
sp.scores.dune.ca <- scores(ca.dune, display="species", choices=c(1,2)) # idem scores axes 1 et 2 seulement
sp.scores.dune.ca <- data.frame(sp.scores.dune.ca)

Comme on a fait acvec les scores des espèces de la PCA, on ne garde que les plus hautes valeurs absolues.

A.ca <- top_n(sp.scores.dune.ca, 3, CA1)
B.ca <- top_n(sp.scores.dune.ca, -3, CA1)
C.ca <- top_n(sp.scores.dune.ca, 3, CA2)
D.ca <- top_n(sp.scores.dune.ca, -3, CA2)
sp.scores.skim.ca <- bind_rows(A.ca,B.ca,C.ca,D.ca)
sp.scores.skim.ca <- distinct(sp.scores.skim.ca)

Graphique simple avec les fonctions plot et ordiplot.

plot(ca.dune)
ordiellipse(ca.dune, dune.env$Management, display = "sites", kind = "sd", label = T)

ggplot

On commence le graphique ggplot. On peut faire un graphique qui ressemble à celui fait par plot en ajoutant “à la main” les sites et les espèces (en fonction de leurs scores extraits).

ca.plot <- 
  ggplot(data = site.scores.ca, aes(x = CA1, y = CA2)) +
  expand_limits(x=c(-4,4), y=c(-1,3)) + # sert à modifier l'échelle numérique du graphique
  geom_text(data = site.scores.dune.ca, aes(x = CA1, y = CA2, label= rownames(site.scores.dune.ca)), color = "black") + # ajout du texte des sites
  geom_text(data = sp.scores.dune.ca, aes(x = CA1, y = CA2, label= rownames(sp.scores.dune.ca)), color = "red") + # ajout du texte des espèces
  xlab("CA1 (25.34%)") + ylab("CA2 (18.92%)") # contrairement avec la fonction ggord, on écrit la variation expliqué par les axes soi même
ca.plot

On peut rajouter les vecteurs de note envfit.

ca.plot <- 
  ggplot(data = site.scores.ca, aes(x = CA1, y = CA2)) +
  expand_limits(x=c(-2,2.5), y=c(-1,3)) + # on zoom in sur l'échelle
  geom_text(data = site.scores.dune.ca, aes(x = CA1, y = CA2, label= rownames(site.scores.dune.ca)), color = "black") +
  geom_text(data = sp.scores.dune.ca, aes(x = CA1, y = CA2, label= rownames(sp.scores.dune.ca)), color = "red") + 
  xlab("CA1 (25.34%)") + ylab("CA2 (18.92%)") +
  geom_segment(data = vect.ca.env, # ajout des vecteurs envfit
               aes(x = 0, xend = CA1, y = 0, yend = CA2),
               arrow = arrow(length = unit(0.75, "cm")), color = "black") +
  geom_text_repel(data = vect.ca.env, aes(x = CA1, y = CA2, label = variable.env, size = 3), color = "firebrick1", fontface = "bold", show.legend = FALSE)
ca.plot

On va mettre des points pour les sites et changer leurs couleurs et leurs symboles en fonction de leur type d’aménagement. On va aussi mettre en graphique uniquement les espèces aux plus grands scores.

ca.plot <- 
  ggplot(data = site.scores.ca, aes(x = CA1, y = CA2), color = Management) + # on ajoute ici arg "color"
  expand_limits(x=c(-2,2.5), y=c(-1,3)) +
  geom_point(aes(shape = Management, color = Management), size=3) + # on change de geom_text à geom_point pour les sites
  xlab("CA1 (25.34%)") + ylab("CA2 (18.92%)") +
  geom_text_repel(data = sp.scores.skim.ca, aes(x = CA1, y = CA2), label=rownames(sp.scores.skim.ca), color = "black") + # on garde geom_text pour les espèces mais on change de dataframe
  geom_segment(data = vect.ca.env, # ajout des vecteurs envfit
               aes(x = 0, xend = CA1, y = 0, yend = CA2),
               arrow = arrow(length = unit(0.75, "cm")), color = "black") +
  geom_text_repel(data = vect.ca.env, aes(x = CA1, y = CA2, label = variable.env, size = 3), color = "firebrick1", fontface = "bold", show.legend = FALSE) +
  scale_shape_manual(values = c(15,16,17,18)) # on décide des symboles à utiliser (se référer à la charte plus haut)
ca.plot

On va ajouter des ellipses autour de nos sites en fonction de leur type d’aménagement. On va d’abord utiliser la fonction ggordiplot de la librairie ggordiplots pour extraire les données des ellipses. (on peut aussi extraire les données de hulls et spiders.) En même temps, ça nous donne un autre beau graphique, mais qui est plus complexe à modifier.

library(ggordiplots)
## Loading required package: glue
ca.ggordiplot <- gg_ordiplot(ca.dune, groups = dune.env$Management, kind = "sd", conf = 0.95, pt.size = 3) # on spécifie l'intervalle de confiance des ellipses

names(ca.ggordiplot)
## [1] "df_ord"      "df_mean.ord" "df_ellipse"  "df_hull"     "df_spiders" 
## [6] "plot"
ca.ellipses <- ca.ggordiplot$df_ellipse # extraction des valeurs des ellipses

On prend les valeurs des ellipses qu’on vient d’extraire pour mettre en graphique dans ggplot ces ellipses. Om change aussi les couleurs de nos groupes.

ca.plot <- 
  ggplot(data = site.scores.ca, aes(x = CA1, y = CA2), color = Management) + 
  expand_limits(x=c(-2,2.5), y=c(-1,3)) +
  geom_point(aes(shape = Management, color = Management), size=3) +
  xlab("CA1 (25.34%)") + ylab("CA2 (18.92%)") +
  geom_text_repel(data = sp.scores.skim.ca, aes(x = CA1, y = CA2), label=rownames(sp.scores.skim.ca), color = "black") +
  geom_segment(data = vect.ca.env, # ajout des vecteurs envfit
               aes(x = 0, xend = CA1, y = 0, yend = CA2),
               arrow = arrow(length = unit(0.75, "cm")), color = "black") +
  geom_text_repel(data = vect.ca.env, aes(x = CA1, y = CA2, label = variable.env, size = 3), color = "firebrick1", fontface = "bold", show.legend = FALSE) +
  scale_shape_manual(values = c(15,16,17,18)) +
  geom_path(data = ca.ellipses, aes(x = x, y = y, color = Group), show.legend = FALSE) + # ajout des ellipses
  scale_color_brewer(palette="Dark2") # changement des couleurs des ellipses et des sites
ca.plot

Enfin, on pourrait rajouter notre thème personnalisé, mais optons pour un existant comme “bw”.

ca.plot + theme_bw()

Comme mentionnée précédemment, vous pouvez modifier presque à l’infini votre graphique!

On sauvegarde notre figure de CA.

ca.plot + theme_bw()

ggsave("ca.png", path="./figures/")
## Saving 7 x 5 in image

PCoA

Préparation des données

Faire une PCoA avec la fonction prcomp.

dune.bray <- vegdist(dune, method='bray')

Faire une PCoA dans R

pcoa.dune <- cmdscale(dune.bray, k =(nrow(dune) - 1), eig = TRUE)
## Warning in cmdscale(dune.bray, k = (nrow(dune) - 1), eig = TRUE): only 14 of the
## first 19 eigenvalues are > 0
summary(pcoa.dune)
##        Length Class  Mode   
## points 280    -none- numeric
## eig     20    -none- numeric
## x        0    -none- NULL   
## ac       1    -none- numeric
## GOF      2    -none- numeric

Extraire les scores de la PCoA.

site.scores.pcoa <- scores(pcoa.dune) # site scores
sp.scores.pcoa <- wascores(pcoa.dune$points, dune) # species scores 
site.scores.pcoa <- data.frame(site.scores.pcoa)
sp.scores.pcoa <- data.frame(sp.scores.pcoa)
Management <- dune.env$Management
site.scores.pcoa <- cbind(site.scores.pcoa, Management)

Comme on a fait acec les scores des espèces de la PCA et de la CA, on ne garde que les plus hautes valeurs absolues, mais cette fois ci, des axes 1 et 3 (pour le plaisir)!.

A.pcoa <- top_n(sp.scores.pcoa, 3, X1)
B.pcoa <- top_n(sp.scores.pcoa, -3, X1)
C.pcoa<- top_n(sp.scores.pcoa, 3, X3)
D.pcoa <- top_n(sp.scores.pcoa, -3, X3)
sp.scores.skim.pcoa <- bind_rows(A.pcoa,B.pcoa,C.pcoa,D.pcoa)
sp.scores.skim.pcoa <- distinct(sp.scores.skim.pcoa)

envfit SUR LES AXES 1 ET 3 (on n’oublie pas!)

(pcoa.envfit <- envfit(pcoa.dune, dune.env, choices=c(1,3)))
## 
## ***VECTORS
## 
##         Dim1      Dim3    r2 Pr(>r)  
## A1 0.9999700 0.0078085 0.379  0.021 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
## 
## ***FACTORS:
## 
## Centroids:
##                 Dim1    Dim3
## Moisture1    -0.2650  0.0607
## Moisture2    -0.1620  0.0000
## Moisture4     0.1065 -0.2369
## Moisture5     0.3271  0.0069
## ManagementBF -0.2694  0.0706
## ManagementHF -0.1350 -0.0394
## ManagementNM  0.1822  0.0370
## ManagementSF  0.0650 -0.0395
## UseHayfield  -0.0608 -0.0240
## UseHaypastu  -0.0227 -0.0210
## UsePasture    0.1213  0.0673
## Manure0       0.1822  0.0370
## Manure1      -0.1654  0.0175
## Manure2      -0.1648 -0.0944
## Manure3       0.1397 -0.0451
## Manure4      -0.1656  0.0944
## 
## Goodness of fit:
##                r2 Pr(>r)    
## Moisture   0.6918  0.001 ***
## Management 0.2634  0.141    
## Use        0.0614  0.640    
## Manure     0.2892  0.203    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999

Pour incorporer les vecteurs résultant de ce test, il faut d’abord les extraire et les mettre dans un dataframe.

vect.pcoa.env <- data.frame(pcoa.envfit[["vectors"]][["arrows"]])
vect.pcoa.env$variable.env <- rownames(vect.pcoa.env)

PERMANOVA. On va tester si la composition végétale diffère entre différents type d’aménagement (facteur; variable quantitative) des sites.

dune.dist.bray <- vegdist(dune, method='bray')
disper.dune.bray <- betadisper(dune.dist.bray, dune.env$Management)
anova(disper.dune.bray)
## Analysis of Variance Table
## 
## Response: Distances
##           Df  Sum Sq  Mean Sq F value Pr(>F)
## Groups     3 0.13831 0.046104  1.9506 0.1622
## Residuals 16 0.37816 0.023635
permanova.dune.pcoa <- adonis2(dune.dist.bray ~ Management, data = dune.env)
permanova.dune.pcoa
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = dune.dist.bray ~ Management, data = dune.env)
##            Df SumOfSqs      R2      F Pr(>F)    
## Management  3   1.4686 0.34161 2.7672  0.001 ***
## Residual   16   2.8304 0.65839                  
## Total      19   4.2990 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Extraction des ellipses de la PCoA (axes 1 et 3!)

library(ggordiplots)
pcoa.ggordiplot <- gg_ordiplot(pcoa.dune, choices = c(1,3), groups = dune.env$Management, kind = "sd", conf = 0.95, pt.size = 3)

names(pcoa.ggordiplot)
## [1] "df_ord"      "df_mean.ord" "df_ellipse"  "df_hull"     "df_spiders" 
## [6] "plot"
pcoa.ellipses <- pcoa.ggordiplot$df_ellipse

ggplot

pcoa.plot <- 
  ggplot(data = site.scores.pcoa, aes(x = Dim1, y = Dim3), color = Management) + 
  expand_limits(x=c(-0.35,0.6), y=c(-0.35,0.35)) +
  geom_point(aes(shape = Management, color = Management), size=3) +
  xlab("PCoA1 (X%)") + ylab("PCoA3 (X%)") +
  geom_text_repel(data = sp.scores.skim.pcoa, aes(x = X1, y = X3), label=rownames(sp.scores.skim.pcoa), color = "black") +
  geom_segment(data = vect.pcoa.env, 
               aes(x = 0, xend = Dim1, y = 0, yend = Dim3),
               arrow = arrow(length = unit(0.75, "cm")), color = "black") +
  geom_text_repel(data = vect.pcoa.env, aes(x = Dim1, y = Dim3, label = variable.env, size = 3), color = "firebrick1", fontface = "bold", show.legend = FALSE) +
  scale_shape_manual(values = c(15,16,17,18)) +
  geom_path(data = pcoa.ellipses, aes(x = x, y = y, color = Group), show.legend = FALSE) + 
  scale_color_brewer(palette="Dark2")
pcoa.plot

ggsave("pcoa.png", path="./figures/")
## Saving 7 x 5 in image
  • FIN -